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The goal of this investigation was to contribute to 
the understanding of natural convection effects in phase 
change thermal control devices. This goal was accomplished 
by developing a mathematical model to evaluate natural con- 
vection effects in a phase change test cell undergoing 
solidification and then evaluating the model against experi- 
mental data. Although natural convection effects would be 
minimized in flight spacecraft, all phase change devices 
would be ground tested and thus understanding the effect 
of natural convection on such devices and the ability to 
predict ground-based system thermal performance become 
quit.e_j.mpor tan t . 

The mathematical approach to the problem was to. first 
develop a transient-two-dimensional conduction heat transfer 
model for the solidification of a normal paraffin of finite 
geometry. Next, a transient two-dimensional model was devel- 
oped for the solidification of the same paraffin by a com- 
bined conduction-natural-convection heat transfer model. 
Throughout the study, n-hoxadecane (n-C^H^) was used as 
the phase-change material in both tne theoretical and the 
experimental work. The models were based 021 the transient 
two-dimensional finite difference solutions of the energy, 

continuity, and momentum equations. The convection model 
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assumed incompressible flow except as modified in the gravity 
terms in the equations of motion. 

An experimental system was set up to verify the 
theoretical analyses and results. The system consisted 
of a closed rectangular box inclined at various angles from 
the horizontal plane, and cooled from below. The box was 
completely filled vrith n-hexadecane , a long- 

chain paraffin. 

Gravity levels were calculated, depending on the angle 
of inclination of the test cell. Temperatures in various 
parts of the cell were recorded by 24 thermocouples as func- 
tions of elapsed time from the start of the experiments. 
Comparisons we re made between experimental results and 
computer-calculated theoretical results, 

Heat transfer when the cell was in the horizontal 
plane was by conduction, Uith_the cell inclined at various 
angles, the heat transfer rate was increased due to com- 
bined conduction-convection heat transfer. The cell was 
generally cooled below the conduction temperatures, when 
convection was also present. The shape of the phase inter- 
face between the liquid and the solid phases was also changed 
l rom a flat plane to a curved surface oy convection. One 
half of the test cell was cooled faster than the other half, 
during convective cooling. Convection was found to be an 
important parameter in the solidification process. 
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Remarkably goou agreement wan obtained between experi- 
mental data for solidification with conduction heat transfer 
and calculated results from the theoretical transient two- 
dimensional conduct ion-hcat-transfor model. Close agreement 
was also obtained between experimental data for solidification 
with the test cell inclined at various angles and the theoret- 
ical results from the transient two-dimensional conduction- 
natural-convection heat transfer model. The trend of the 
effect of "convection on the system was clearly revealed by 
both the experimental -and the theoretical results. 

Stability problems ’were encountered in the finite- 
difference solutions of the theoretical heat- transfer models. 
Central differences were used to eliminate the dependence of 
the finite-difference time step on the unknown velocity com- 
ponents in the convection equations which would have occurred 
if forward or backward differences had been used. Maximum 
allowable velocities were severely restricted by stability 
requirements of the finite differences , so that we could not 
use as high velocities as we would have, liked to in the solu- 
tions of the convection-temperature equations. In addition 
to these restrictions , more severe restrictions on the magni- 
tude of a time step were imposed by the need to minimize 
errors due to numerical dispersion on the theoretical results _ 
from the natural convection model. Numerical dispersion error 
terms, introduced by the neglect of second-order time deriv- 
atives in the central difference approximations of the bound- 
ary layer energy equations, could not be eliminated entirely, 
but they were minimized to less than 10 percent by the use of 
very small tine-steps. 

This study is a first attempt to model the effects of 
natural convection on a solidifying system with moving bound- 
aries. It is also a good basis for further, more rigorous 
analyses in which the energy, momentum and continuity equa- 
tions may be solved completely. 
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INTRODUCTION 


Phase-change phenomena have received wide scientific 
attention for some time and are of significant importance 
in many technical problems such as solidification of a 
billet, formation of snow, solidification of an asphalt 
layer, melting of alloys and growth of crystals. Recently, 
phase-change materials have been seriously considered for 
spacecraft thermal control. In concept, such materials 
would be used in passive systems employing the process of 
melting or solidification to remove or add themal energy 
from or to a system. Thermal control systems based on 
solid-liquid phase change have many advantages which make 
them useful for certain applications. They are light, easy 
to handle, and easily used as wall-lining elements around 
electronic equipment. 

When a phase-change thermal control material goes 
from the solid to the liquid phase or vice versa, convection 
currents may be set up in the liquid phase by temperature 
gradients in the system. These temperature gradients affect 
the density of the material and hence the force of gravity 
on it. Thus, the resultant convection is known as gravity- 
induced or natural convection. It is not a forced convec- 
tion since no initial bulk flow is forced on the system. 
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The present investigation was aimed at studying the 
effect and importance of such natural convection currents 
on the rate of solidification, temperature profiles and 
the shape of the interface of a phase-change material, 
n-he?:adecane (n-C^H-^) , a normal long-chain paraffin that _ 
is liquid at normal room temperatures* It was hoped that 
the investigation would help reveal any influence of natural 
convection on the ability of such a phase -change material 
to control the temperature of an instrument around which 
it is packaged* If convection increased heat transfer 
remarkably, a_decision would have to be made as to -the 
suitability of using such a material in_the liquid phase 
for thermal control. 


I-.uch theoretical work has been presented in the 
literature dealing uith problems which are related to 
physical change of state. The basic feature of such pro- 
blems is the existence of a moving boundary or surface 
between phases. Therefore , the problem that is most often 
considered is how to determine the way in which this sur- 
face or boundary moves. Heat may be liberated_or_ab sorbed 
on the boundary ; there may be a volume change accompanying 
the change of state, and the thermal properties of the pha.se 
on either side of the interface may be different for the 
phases and may vary as the change of state proceeds. There- 
fore, the problem is non-linear in nature and general ana- 
lytical solutions for it are not available. Some exact 
solution* for models that mathematically approximate the 
real problems have been obtained, mostly for infinite or 
semi -inf ini te geometry. 

Car slav; and Jaeger ^ were among the first to give 
in-depth treatment of melting and solidification problems. 
They commented on the need to use numerical and finite- 
difference methods in solving many of the complex problems 
that arise in finite geometric configurations. 

Many of the solutions presented in the literature 
of phase-change problems are valid only if the material 
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under study is initially at. i ta couilj.br ium temperature 
for chcuigc of state* — They ignore the most frequently 
encountered case in which the. material under study is 
initially at a temperature quite different from the phase- 

change temperature. 

( 2 ) 

Stefan was the first to give a published discussion 
of a one -dimensional transient conduction problem vdth phase 
change, for a single component or eutectic composition vdth 
constant properties. Thus, the term "Stefan's Problem" 
came to be used to describe a one -dimensional conduction _ 
problem in which a semi-infinite slab initially at a constant 
non-zero temperature, has one face maintained at zero temper- 
ature for time greater than zero. The solution to the pro- 
blem entailed the assumption that the time-dependent interface 
position was proportional to the square root of the product 
of time arid the thermal diffusivity cf the material of the 
slab. 

Danckwerts^ , Booth^^, and Kreith and Bomie^^ have 
all presented analytical or semi-analytical solutions to 
phase-change problems under various boundary and initial 
conditions. Chao and Weiner^ investigated the temperature 
in a solid-liquid system while the liquid was being poured. 

The latent heat of phase change was treated as a "pseudo" - 
specific heat and the solution, obtained by a Laplace trans- 
form technique was an integral that v/as solved numerically. 
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Kany authors used tho variational technique to solve- 
hcat transfer problems, with or without a phase change. 
Chambers Biot and Daughadayi^ used this approach. 

Biot and Daughaday studied an ablation problem in which 
the melt was removed as it was formed. 

The heat-balance-integral technique, an analytical 
method that gives approximate solutions to a wide variety 
of heat transfer problems, is used in many papers, in the 
literature. It is mostly used to solve non-linear prpblems 
that must be solved numerically or approximately. Its big 
advantage is that it changes an energy equation from a partial 
differential equation to an ordinary differential equation. 

One disadvantage of this method is that the solutions ob- 
tained satisfy the differential equations only on the average. 
Goodman^) and Poots^ 10 ^ have used this method to study heat 
transfer problems. — Poots studied the two-dimensional inward 
solidification of a uniform prism initially at the fusion 
temperature. 

In the study of more general cases of phase-change 
problems, numerical analysis may be the only feasible tech- 
nique available. Dusinberre^ 11 ^ , Pujado^ 12 ^, Ukanwa, 

S ter mole and Golden (1 ^^ have all used finite-difference 
techniques to study phase-change problems in which the heat 

transfer mode is by conduction only. 

( 14 ) 

Killer' used the "surplus temperature" technique 
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in an attempt to improve the predictions of the phase front. 
To account for the heat absorbed at the phase front, the 
calculated temperature was permitted to exceed the actual 
melting temperature until an arbitrarily selected tempera-_ 
ture was reached. Whan this temperature was reached, the 
grid element containing this particular, nodal point was 
considered to have changed phase, and the phase front was 
shifted to the next node. 

Ehrlich^. 1 gave the implicit finite-difference 
equations for the one-dimensional melting problem with a 
variable heat input specified as a function of time. The 
implicit equations were then put into tridiagonal matrix 
forms for easy solution by Gauss elimination and by back- 
substitution. Special modified equations were given for 
nodes near the freezing front. 

The Northrop Corporation reports^ 1 ®' ^ presented 
a survey of the phase-change problems involving selection 
of the proper compounds for phase-change thermal control 
devices, evaluation of properties, and experimental study 
of different test cells. Some of the physical properties 
used in the present study on n-hexadecane have been taken 
from these reports* 

Other works on phase-change phenomena include Bannister 
and Bentilla^ 1 ®) , Grodzka and Fan^ 1 ^, and Chambre^ 20 ^. A 
survey of many papers on phase-change phenomena has been 
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presented by Huehlbaucr and Sunderland. 

Convection in enclosed fluids has been studied exten- - 
sively. Unfortunately, most of the studies have been either 
only experimental with no attempt at theoretical modeling 
or only theoretical with no experimental corroboration. In 
addition, the majority of the studies on convection deal 
with flat or parallel plates and on gases. The works that 
do deal with completely- enclosed liquids often have theo- 
retical solutions that hardly, if ever, approximate real 
situations closely. 

Some good texts for theoretical references on heat trans- 
fer and fluid flow are Rohsenow and Choi' 22) , Schllchting (23) , 
Longwell' 2 ^), and Bird, Stewart and Lightfoot' 25 ). The last 
two references- were used as the .references for developing 
the basic boundary-layer equations for gravity^induced oon- 
vection in the present study. 

Uilkes and Churchill' 26 ) studied temperature profiles 
in an enclosed rectangular cavity subject to adverse temper- 
ature gradients. Their equations for gravity- induced free 
convection were developed from the basic equations of con- 
tinuity, motion and energy. The resultant system of equations 
was solved by an implicit alternating-direction technique 
developed by Peaceman and Rachford.^^ 

Chandrasekhar' 28 > gave an extensive treatment of 
stability and instability in fluids subject to adverse 
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temperature, gravity and magnetic effects, .'l’he method used 
by Pellew and Southwell linearize temperature and 
velocity equations was.-also discussed. It involved -the 
solution by the separation-of-variables technique, the 
equations - of continuity, momentum and energy in order to 
determine the critical Rayleigh number (gh-^dd A T ud T ) xiec- 
essary to initiate free convection. in a fluid heated from 
oeiow. The method assumed that motion in fluids heated from 
below was cellular and involved finite numoers of rolls 
corresponding to particular wave numbers. 

Leont'ev ana Kirdyasmsin^-^ studied convection m 
fluids of large volumes. Tney assumed tnat, except in tne 
boundary layer which was very small relative to the dimen- 
sions of the fluid, motion was by ideal flow. However, the 
maximum velocity, occurred in the boundary layer, near the 
walls, where the temperature gradient or heat flux was 
largest. This ideal-flow approx ima t ion_£or natural convec- 
tion in completely enclosed fluids has been used by other 
authors. The velocity profiles used in the present study 
assumed ideal flow profiles similar to their models. This 
approach was necessitated by the difficulty in satisfying 
all no-slip conditions on rigid boundaries of completely 

enclosed fluids. 

( 31 ) 

Bain w performed a two-dimensional experimental 
study of gravity-induced convection in an enclosed liquid 
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in a rectangular cavity inclined at an angle from the !tior- 
izontal plane. Solid-liquid phase, change was also involved. 
Person a l communication with Mr. Bain has revealed that. he 
is presently involved in providing theoretical ideal-flow 
models to predict his experimental data. 


THEORETICAL A^JALYS IS 


Formulation of the Problem 

The_physical problem to be studied is the solidifica- 
tion of normal hexadecane enclosed in a cell cavity of 
height h (initial height h Q ) and constant square cross- 
sectional area of lateral dimension W. The cell is inclined 
at angle a degrees measured from the horizontal plane. The 
coordinate axes and their origin are located as shorn in 
Figure 1 with x along the cold bottom plate and y perpen- 
dicular to -it. 

Some physical properties of the test material, n- 
hexadecane (n-C^H^) , were taken from Reference 17: 

Density 

Solid n-hexadecane: d s = 1.0772-8.41 x 10" 4 T gm/cm 3 
for T L 290. 0°K “ 

Liquid n-hexadecane: d L =0.9726-6.8l3 x 10" 4 T gm/cm 3 
for 290. 0°K L T L 400. 0°K 

Specific Heat 

Solid n-hexadeane: c pS =0.50 5 cal/(gm-°K ) 
for T L 290. 0°K 

Liquid n-hexadecane: c pL =0. 1626+1 .164 x 10“ 3 T cal/(gm-°K) 
for 290. 0°K L t L 480. 0°IC 
Thermal Conductivity 

Solid: k s =2.390 x i0" 3 -3.047 x 10" 6 T watt/(cra-°K) 
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Liquid; 1< L =2.390 x 10” 3 - 3.04? x 10“ 6 T watt/(cm-°K) 
for 250°K L T L 425°K 

Solidification temperature; Tp = 290. 0°K = l6.7°C 
Latent heat of solidification; H f = 56.67 cal/gm 

The cold bottom plate is maintained at temperature, 
T(x,0, t) = f(t). The effects of natural convection,, induced 
by density and temperature changes, on the solidification 
rate and the temperature profiles of the paraffin are to 
be studied as functions of time and angle (equivalent to 
changing gravity levels). By changing angle a, one may 
vary the components of the force of gravity acting on the 
system. Only density changes affecting the gravity term in 
the equations of motion are considered. 

For this investigation, temperature-averaged properties 
are used. The following properties, obtained by experiment 
and from References 17 , 12 and 33 are used: 

Tf = 290.66°K (experiment) 

dg = O.833 gm/ciP (experiment at 268. 2°K) 

85 0.755 gm/cm^ (experiment at 294. 7°K) 
dp = density of plexiglas (material of cell wall) 

55 1*155 sm/cm^ (experiment at 298°K) 

From Reference 1?: 

3 = -6.813 x 10"^ gm/(cm^-°K) (from d^=d + jJT) 
c pL = 0.506 cal/(gm-°K) (at 296. 0°K, the average temp- 
erature encountered in the liquid phase) 
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Opg- 0.505 cal/(gm-°K) 

k S = 3.720 x 10* 4 cal/ ( sec-cm-°K) (at 280. 2°K, the 
average temperature encountered In the solid phase) 

k L = 3.555 x 10* 4 cal/(sec-cm-°K) (at 296 . 0°K) 

Prom Reference 32: 

k p - plexiglas thermal conductivity = 4.960 x 10* 4 
cal/(sec-cm-°K) 

c pP ss Pl GX iSlas specific heat = O.35 cal/(gm-°K) 

Prom Reference 33_: 

v = kinematic viscosity of n-hexadecane =0.04 cm 2 /sec 

Prediction of the rate of solidification and the 
temperature profile of the test paraffin under a given set 
of initial and boundary conditions are to be made using 
(l) a conduction model and (23 a combined conduction-convec- 
tion heat-transfer model. The following definitions are 
used? For a function f(x,y,t), 

6f/6t = f t ; 6 2 f/6t 2 = f t t ; a f/Sx = f x ; jZf/** 2 -^! 


6f /*y = f y ! « 2 f/«y 2 « fyy ; « 2 f/6x6y - f xy (!) 

Df/Dt = f t + Uf x + vfy (2) 

With these definitions, we have 

(1) Conduction Model ; temperature equations: 

Liquid phase 

T t = X L (T xx + T yy } f0r Y(x » t) ^ L * L * 0 (3) 

T(x,Yj t) = T f (i) 
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T (*»h 0 + h p ,t) = T a (11) 

Jy(o.y.t) " ^O'.y.t) = 0 dll) 

Kx,y,0) = T a (iv) 

Solid phase 

T t = X S (T xx + T yy> for 0L y L Ux,t) (if) 

T(x,Y,t) = T f (i) 

T(x,0, t) « f(t) (ii) 

T x (0,y,t) s T x (W,y,t) » 0 (iii) 

T(x,y,0) = T a ( iv ) 


In Equations (3) and (4), T a is ambient temperature, h 
is the thickness of the plexiglas wall of the test cell 
and X = k/(dCp) is the thermal diffusivity of the test 

material, ho equations of motion are involved in the con- 
duction model. 

(2) Combined Conduction-Convection Model 
Temper at ure equations: 

Solid phase:- Equations (4) apply. 

Liquid phase 

DT/Dt = (T xx + T yy > ( 5 ) 

Viscous heating is negligible. The boundary conditions (i) to 
(iv) of Equation (3) apply here also. 

Equation of Continuity: 


( 6 ) 
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for mi incompressible fluid. 

Equations of Motion: _ 

The Navier-Stokes' equations of motion for an incompressi- 
ble fluid are 

d L Du/Dt = d L i> ( Uxx + Uyy) - P x -_d L g sin o (?a) 

d L Dv/Dt =■ d L v (t xx + Vyy ) - P y - d L g cos a (?b). . 

Velocity components u and v are each zero at each of the 
walls, x = 0, x = V/, y ■ 0, and y * h Q .. At time t = 0, 
velocities are also equal to zero. 

The temperature equations for both liquid and so- 
lid phases are coupled by mi interface equation for phase 
change involving enthalpy change during solidification. This 
is true for both models: conduction and combined conduction- 
convection models. The proper phase-change equation will be 
discussed later. 

Equations (7) may be modified as follows. Let 
P - P + p, where P is the value of the pressure in the test 
cell when there is no motion. Also, T and d L are the values 
of temperature and the density of the test material (liquid) 
when there is no motion. V.'hen there is no motion, Equations 
(7) become 

0 * - 7 X - sin a (7»a) 

0 * - Py - d L g cos a (7»b) 

Uhen Equation (7* a) is subtracted from Equation (7a),. and 
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(7'b) from (7b), the following modified equations are obt- 
ained for gravity-induced natural convection as a result of 
temperature-caused density variations in the test material. 

Du/Dt « v (u^ + Uy y ) - p x /d L - pg(T-T) sin a (8a) 

\ 

Dv/Dt = v ( Vxx + v yy ) - Py/'d L - pjs(T-T) cos a (8b) 


where 3(T-T) = d L - d L » Equations (.0) inolude_the effect of 
the angle of inclination of the test cell from the horizontal 
plane. Thus the normal gravity level is modified by the angle, 


Equations (8) may be further modified as follows! 
Let a stream function, </>(x, y,t>, be defined by the two 
equations 


u = “ 
v = © 


(9a) 

(9b) 


Thus, Equations (9) automatically satisfy the continuity 
equation (6). Define a vorticity. 6i>(x.y f t) f by the two 
identical equations 


" s v x - u y (10a) 

03 ~ ^xx + ^yy (10b) 

Differentiate Equation (8b) with respect to x and Equation 
(8a) with respect to y, and subtract the result of the later 
operation from the result of the former. The resultant 
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equation becomes , on tho ap plication of Equation (6) end 
Equation (10); 

Do) /Dt * r i fJ L x + ^ yy ) + if ( (T~T) y sin a - (T-T) x cos a) (11) 

d L .. 

Thus, Equation (11) r eplae e s~Equat Ions. ($) as the equation 
of motion. The pressure terms drop out of the equations. 

As is shown in Appendix C, limitations imposed by . 
stability criteria for stable solutions by finite-difference 
approximations of Equations (5) and (11) or (5) and (8) 
simultaneously could not be satisfied within the limitations 
of digital -computer time and memory available to us. In 
order to get meaningful results from these full equations. 
of motion and-Jtemperature, one would have to use very small, 
time and space grid elements, values that are too small to 
allow for the acquisition of meaningful theoretical data 
within. the time and memory limitations permitted to these 
investigators by the university computer center. _ 

Therefore, tno approach to be used to obtain approx- 
imate solutions of Equations (5) and (8) is to use approxi- 
mate velocity profiles obtained by combining an ideal-flow 
model with a viscous-flow model. 

Approximate Voice 1 tv Profiles 

Assume that (T-T) x is negligible, compared to . 

y 

(a) Ideal-flow Model (Maintained Convection) 
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Let T-T « c u(x,y,t) where o is a constant, and let 
= '“ <l ^(x^t) where q is a characteristic constant. 
Also, assume that </>(::, y,t) is separable into a product of 
functions of x, y and t. Prom Equation (JL1) , W e gat 

( 12 ) 


Vt - + V + HE ty sin a 

qd L 


ivith initial and-boundary conditions,. 

^(0»yit) = <p (W,y , t) =_<p(x,0,t) 33 -</>(jx.,h, t) . (13) 

V/hen Equation _( 12) is solvcd_by the separation-of -variables 
method, we get 

^n,m =A n,m 5lll( ~ )T,X e -n( q *f (2m+D 2 ggQ| 

^ qd T rh2 


2 ( 2n+l ) ^ (2m+l) 

where q = rr^( — -- + — 

VT h 2 


2 . _ 


2 ), m 0,1,2, • . . ; n=0 ,1,2,.., 


(141 


For maintained motion (independent of time), the exponent 

in the exponential term is zero, that is 

c = -d L v q 2 h 2 

T" T~“ (15) 

ir 4 (2m+l)*0g 

When Equation (15) is put into Equation (14), we get 

’n.m = V ra 3ln 3ln ( 2m+l)w y 

V/ h 1 

mss 0,l,2, ..... ; n=0,l,2, 

Equation (16) is analogous to the equation for the vibrations 
of a striae in which q (Equation 14) G ives the character- 
istic numbers correspondins to the modes. If n = m = 0, we 
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Cot a single circulation flow (one finite roll) in the test 
cell. It either n or m or both are not equal to zero, mult- 
iple rolls-rosult* Figure 2 chows sketches of the rolls that 
would result for various values of m and n, 

(b) Viscous-flow liodel (First-order Kantorovich profile 
approximation by variational method). 

Assume that T-T = yiT Q /h , where bT Q is the temp- 
erature difference betvreen the hot plate at y = h and the 
cold plate at y = 0. ALse, assume a steady-state solution 
and negligible pressure variation in Equation (8a). There- 
fore, Equation (8a), as modified, becomes _ 


u + 
XX 


u, 


yy 


i3gyuT- 

— sin a 


0 


d L «>h 

Assume a profile , (pm h(y/h - 2y 3 /h 3 + yVh 4 ) X(x) 


(17) 

(18) _ 


where the profile has been made to satisfy in the y-direction 
the boundary conditions of Equation (12) for a stream func- 
tion. From Equations (9a) and (18), we get 

u - (1 -6y 2 / h 2 + Ay 3 /h 3 ) X(x) ( 19 ) 


The function X(x) is to be determined by variational technique 
and should satisfy the boundary conditions 


X(0) = X(W) * 0 

The appropriate equation for X(x) 


( 20 ) 

by variational method is 



\ 


obtained by integrating the product of l-6y 2 /h 2 +4y3/h 3 and 
Equation (17) with respect to y between y = 0 and y =_h. 
The definition of Equation ( 19 ) for u is employed in the 
integration. The resulting equation for X(xL is 


d 2 X/ dx 2 - (168/I7h 2 ) X + . 

— ■ — ■ sin a = 0 
17 d T 


( 21 ) 


On solving Equation (21) with the boundary conditions of 
Equation ( 20 ), we get 

y - £gh 2 -iT. sin a . . ... . ± , 

x 0 /, _ sinh y(U~x)+sinh yx \ /ooX 

2^i’d — — ~ 


sinh yVl 


where y = Vt 1^37177 /h 

Therefor, the first-order Kantorovich profiles are 


(23) 


u = -B(l-6y 2 /h 2 +4y 3 /h 3 )(:i- s lDh V(W-*>*slnh yx, (24a) 

sinh yW 

v = Byh(y/H-2y 3 /h W) i—jS Ux J-° osh YX ) (24b) 

' sinh yW / 

tn -nut /u r> 3 /, 3 i± „ A. , Y(W-x)+sinh yx 

9 - Bh(y/h-2yyh 3 +y yh^) (1- rya-m ) (24c) 


whore B = -|!gh 2 CT 0 sin a/(24i’d, ) . 


The stream function of Equation (24c) gives one finite roll 
(singlet flow) * As it turned out, however, a single symme- 
tric roll could not predict the experimental data in all 
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sections of the tostjcelL. It could predict only one half 
of the test cell (QdxCW/2) • Multiple rolls given- by Equa- 
tion (1 6) when at least one of m and n is not zero, do not 
predict the entire cell performance either. 

Hence, a. modified profile involving a combination 

of Equation (24c) and Equation (16) for. n = m * 0 was used 
to get a single roll, non-symmetric in the-x-direction, but 
symmetric in the y-direction. . Equation_(l6) , for n = m * 0, 
was allowed to apply for Q 6 x L VP/2, so that at x = VP/2, 

v = 0; and u = --A 0|0 g cos gy (25) 

Equation (25) was obtained by applying Equation (9a) on 
Equation (16) and letting x = Vi/2. VP/2 must satisfy the 
relationship 2h 1 Vi /£ 1 Vi* The second portion of the 
stream function profile was obtained by using Equation (24c) 
in the. range VP/2 £ x L Vi. In addition, v from this profile 
was set equal to zero at x ~ VP/2 so that u for this part, 
became 

u « -D(l-6y 2 /h 2 +4y 3 /h^) at x=W'/2 (26) 

Equations (25) and (26) should yield the same values for 
u at x=’P/2. They both go to a maximum, a minimum and zero 
at the same points: y/h * 1, 0, £, respectively. Thus, the 
maximum value given by Equation (25) was set equal to the 


f 
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maximum value given by Equation (26) to get A Q c *0h/iT. 
Thus, the complete stream function that was used was 


Bh sin irx sin 1Ty 

i "17 1 h" * for 0 ^ x ^ W»/2 

9 = Eh(yAl-2y 3 /h 3 +y 4 /h 4 ! (l- sinh Y<W-r)+sinh y(W+x-W), 

sin.hr v"C2u:un 

for W«/2 A x A w (2?) 


Figure 3 shows graphs of u from Equation (25) and u from. 
Equation (26) at the merging of the two portions of the 
stream function in Equation (27), that is, at x = w*/2. 

For the present study, W*/2 was chosen as* 0.95 W. 
So, the final form of the resultant stream function profile 
that was used in this study is 


— sin ux si nW 

ir I^w 5T* for 0 t x k 0.95 W 

(p = 

Bh(y/h-2y 3 /h 3 +yVh 4 ) (l- s inh Y (W-x)+sinh yU-0.9 W)^ 

sinh 0.1 yw 

for 0.95 W i x £ V/ (28) 

The introduction of the definitions in Equation (9) into 
Equation (28) gives u and v. Hence, Equation (28) is used 
in the finite-difference solution of Equation (5) to obtain 
the convection temperature profiles in the liquid phase. 
Since the angle, a t appears only in the term B, Equation 
(28) indicates that the angle of inclination, a, affects 








1 


25 


the magnitude of the velocities but not their shape. In 
Equation (28)_, h is the time-dependent height of the liquid 
phase which becomes smaller than the height h Q of the test 
cell cavity as solidification continues. 


i I f orcnc e Approximations of Governing Temperature 
Eauations 


The following definitions are used. Ignore third- 
order derivatives. Figure 4 explains the space— grid arrange- 


ment. 

Central**diff erence approximations 
T ( i+1 , j ) - T(i-l,j) 


2 Ax 

T( 1,3+1.) - T(iJ-l) 

TSy 



(291) 

(29ii) 


T(i+l,J) - 2T(i, j) + T(i-1 , j ) 

(Ax ) 2 * T xx 

T(i, j+1) - 2T(i,j) + T(i, j-1 ) 

— — 5 — = t 

(Ay) 2 yy 


(29iii) 

(29iv) 


Forward-difference approximations 
T'(i.j) - T(i,j) 


at 

2(1+1, J> - T(i| j) 


= T t + 4^ti tt 


(29V) 


ax 


= T x + iaxT^ 


(29vi) 
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T( 1,3+D - T(l, J) 


Ay 


^ T y + 


Backward-difference appi 


(29vii) 


=Ix .^ Txx 

(22.viii) 

- g( i±ii - T y - iAy Tyy 
Ay 

(29ix) 

Coefficients: 


G = 1 - 2XAt(l/(Ax) 2 + l/(Ay) 2 ) 
1 

(30i) 

G 2 = XAt/(Ax) 2 

(30JUJL- 

g 3 = XAt/(Ay) 2 

(30iii) 

C 1 - 2d P C pP h p + d L C pL^ y - - 

V30iv) 

c 2 = kjAy + 2kph p 

C30v) 

G 4 * 1 ^(C 2 /(ax) 2 + k L /oy + kp/h p ) 

(3pvi) 

G^ = AtCg/lC^Ax) 2 ) 

(30vii) 

G 6 * 2Atk L /C t Ay 

(30viii) 

G ? = 2Atk p /C 1 h p 

(30ix) 

Gg * G 2 + u(i,j)At/2Ax 

(311) 

G 9 = G 2 - u(i,3)At/2Ax 

(3iii) 

G 10 * G 3 + v(i,3)^/2Ay 

(31111) 

G 11 * G 3 " v(i,3)At/2Ay 

(3Hv) 
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On the application of tho central-difference, space 

definitions and the time-forward-difference definitions to - - 

Equations (3), (4) and (5)» we get the following sets of 
equations. Second-order time derivative is ignored. 

(1) Conduction Model 
Liquid phase 

T»(i,j) = G iL T(i,J) + G 2L (T(1+1,J) + T(i-lJ)) + 

G 3L (T(i, j+1) + T(i,j-l)),_for 2 £ i £ M-l 
and N3(i)+1 £ 3 £ N-l (32i) 

T'(l,3) = G 1L T(l,j) + 2G 2L T(2J) +..G 3L (T(l,j+l) 

+ T(l,j-1)), for NS(1)+1 L j L N-l (32ii) 
T'(MJ) = G 1L T(MJ) + 2G 2L T(M-1, j) + G 3L (T(M, j+1) 

+— T(M , j-1 ) ) , NS (M) +1 £ 3 L N-l ( 32iii) 

T«(i,N) = G 4 T(i,N) + G^(.T(i+l ,N) + T(i-1,N))+ 

G 7 T a +G 6 T(i » N " 1} » for 2 £ i £ M-l (32iv) 

T«(1,N) = G 4 T(1,N) + 2G 5 T(2,N) + G 6 Tj(1,N-1) 

+ G ? T a (32v) 

T*(M,N) = G^T(M,N) + 2G 5 T(M-1,N) + G 6 T(M,N-1) 

+ G ? T a (32vi) 

Solid phase 

T*(i#5) = G ls T(i,j) + G 2S (T(i+lJ) + T(i-1,3)) 

+ G 3S (T(i,j+i) + T(i,j-D), for 
2 £ i £ M-l, 2 £ j £ NS(i) 


(33D 
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- G 1S T(1,J) + 2Q 2s T(2,j) + 0 3s (T(l,jtl) 

+ T ( 1 , j-1 ) for 2 L j l NS ( 1 ) (33ii) 

= G 1S T(K , j) + ( M —1 , j ) + G^g(T(K,j+l) 

+ for 2 L $ l NS(M) (32iii) 

T 1 ( i 1 1 ) — f ( t ) , for 1 l. i £• M ( 33i v) 

(2) Combined Conduct ion-Convection Model 
Solid phase:- The solid-phase Squat ions (331) to (33iv). of 
the conduction model apply here in their entirety. 

Liquid phase 

T»(i,j) ■ G 1L T(iJ) + G 8L T(i-i,j) 4 G 9L T(i+l,j) + 

G 10L T(i »J" 1) + Gu L T(iJ+l), for 
2 L i L M-l , NS(i)+l L j l n-1 ( 34!) 

Equations (32ii) to (32vi) also apply here in their entirety. 

In Equation (32) to Equation (3 /*) , subscripts L 
and S on a coeificint imply that the_properties of liquid 
and solid test material are to be used respectively, 

Peculat i on of the _Phage^-chan,'Te Rate and the Height of the 
Solid Phase 

In Equations (32) to ( ) , NS(i) is the number of 
y-nodes in the solid phase given as a function of horizon- 
tal location x (that is, tox) . Thus, the height of the 
solid phase at any time t is given by the equation 


r(l,t) a ay(ITS(i)-0.5), for and NS(i)^2 


(35) 
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If N.j(i) ~ 1 , then Y ( 1 , t ) is on the cold bottom plate and is 
equal to Ay/ 2. N3(i) and Y(l,t) must satisfy, at any time, 
tho following inequations. 

1 L NS(i) l N, f or 1 L i c h _( 36i) 

0 ± Y(i,t) £ h, for It i l m ( 36ii) 

At any time t, NS(i) is evaluated as follows. Firstly, 
T(M,tL= f(t) is checked to see if f(t) L T f , The first 
time that f(t) l t^, NS(i) is set equal to 1 for all 
1 L i to see if T(i,NS(i) old + 1) L T f ._ If it is, then a 

new FS(i) is obtained from the equation NS(i) « NS(i) + i 

UJ*vL 

for- that particular i, and thenceforth, that node is 
treated as a node in the solid phase. 

The latent heat of solidification is used to modify 

the specific heat of the solid phase so as to incorporate 

the effect of change of phase. It is assumed that, since 

the rate of change of phase is so slow and the solid phase 

is at a much lower temperature than the liquid phase, most 

of the heat liberated by solidification goes to warm up the 

solid phase in the form of sensible heat. Thus, an effective 

heat capacity Cpg , based on a lumped average temperature, 

T av * fe‘(T.p+T(X|0 , t) ) = a(‘ff+f(t)), of the solid phase is 

defined to include ti*e phase-change enthalpy change; 

Enthalpy change s Latent heat + Sensible heat gain 
per gram per gram per gram 


(37i) 
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Thus, o^(T r 'i? av ) = H f + C pS (T r T av ). (37ii) 

So, we obtain an equation for the effective specific 


heat c 


pS •" 

°ps ” c p S + 2H f /(T r f(t)) • 


(38) 


c pS» instead of c pI ;j, is used in the calculations of the 
temperature profiles of the solid phase in Equation- (30) 
and Equation (33) • No separate interface equation is now 
needed* 

Therefore, to summarize, when the temperature of a 
node in the liquid phase equals- or -drops below the solid- 
ification temperature, T f , the node is henceforth treated 
as a node in the solid phase, and c^ s , instead of c pS , is 
used to calculate solid-phase temperatures so long as 
solidification is in progress. Thus, Equation (28) to 
Equation (38) are sufficient to determine the temperature 
profile, the rate of solidification, and the streamlines 
in the test cell at any time t. 

Stability Criteria for the Finite-Difference Approximations 
(1) Conduction Model 

The stability criteria for this model are 

£ 0 i % k 0 (39) 

for both liquid and solid phases. From Equation (39), maximum 
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alio v/ablc At for given Ax and Ay. are evaluated. Which- 
ever At from Equation (39) is smaller is the controlling - 
^t for a stable solution. It is recommended that. a smaller 
value of At than the maximum allov/able At be used. The 
smaller the value of At, the smaller the error in the 
finite-difference approximation gets. For our system, 
controls the allov/able At. 

(2) Conduction-Convection Model 

The stability criteria of Equation (39) also apply 
here. Other criteria tire imposed on velocity, thus: 

°3 - 0; G 9 - ° ; G 10 > 0; and G n > 0 (40) 

These inequalities yield the following restrictions on u 
and v: 

lu ( i , j )J < 2 X l /Ax , |v(i,J)| £ 2\ L /Ay (41) 

Thus, for a given Ax or Ay, a maximum allowable absolute 
value of u or v may be determined. The resultant velocities 
evaluated for our test material are small, a fact that 
caused tho difficulty in finding any stable finite-differ- 
ence solutions for the full-blown coupled equations of 
energy and motion so long as the gravity term remained in 
the velocity equations. 

The velocities obtained by differentiating Equation 
(28) with respect to x or y must satisfy Equation (41). 
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With the experimentally-observed temperature difference 
between the hot and the cold plates, we could not use the 
full magnitude of Equation (28) without violatingi.Equations 
(41). So, vie restricted the maximum velocities that we used 
to be 95 percent of the maximum velocities pe rmit ted by 
Equations (41), and the constant B in Equation (28) was 
modified accordingly. If the maximum permitted velocities, 
in stead of some slightly lower values were used, our 
solution would be in the critical region. where a small 
round-off error from the digital computer could_make the 
solution unstable and meaningless. Table 1 gives the time 
elements and velocities permitted by stability restrictions 
and the actual values used in our calculations. 


Table 1: 


TM, IIlCHiv^W TS AITD VELOCITIE S ALLOWABLE 'BY STA- 
BILITY CRITERIA 


£x « Ay; \ L = 9.306 x 10 cm 2 /sec 


Space Element 
Ax (cm) 

0.254 

0.127 

OJl£35 


Time Element 
At-(sec) 

17.33 

4.33 

1.08 


u(max) 

(em/sec) 

7.32 x 10"? 
14.65 x 10"? 
29.31 x 10“° 


v(max) 

(cm/sec)L 

2.75 x 10"? 
5.49 x 10"? 
10.99 x 10" J 


Actual values used: 

Ax * 0.127 cm; At = 1.0. sec; u(max) * 13.92 x 10“^ cm/sec; 
v(max) = 5i22 x 10"^ cm/sec 


^ror An alysis of the Centra l Difference Method Used 
(1) Conduction Model 
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The terms ignored in the finite-difference approxi- 
mations are of the order of (Ax) 2 ? , /6, (Ay) 2 T /6, and 

xx yy 1 

i^tT tt . These terms ere assumed to be small provided Ax, 

Ay and At are small. We ran our digital -computer programs 
f or_a square-grid system (Ax — Ay) f or - different values of 
Ax and At, l/e found that a solution at Ax = o T i ? y cm, At = 

2.0 sec was not much different from a solution at Ax * 

O.O635 cm, At_= 1.0 sec (see Table 2), but involved quite 
^ difference in computation time. Hence, we assumed 

that the- solution had converged at Ax = 0.137 cm and At - 

2.0 sec, and therefore, used ^t = 1.0 sec and Ax = 0.127 cm 
for computation. 


Table 2: 


^CILFOB-COiT^p.GIT^ S OP FIKira -DIPPEKEKCP. SOLUTION 
OM THE GOHDUCTlOj MODEL 


Time 

(1000 sec) 


Ax = 0.127 cm 
^t = 2.0 sec 


o ^ 

°K 

Ax = 0.0635 cm 
A t = 1.0 sec 


« O.36 

1.44 


2,16 
2.52 
2.83 
3.60 
• 96 
.63 
5.04 
5.76 
6.48 
7.20 


l 


299.636 
293.139 
291.272 
290.622 
290.185 
289.472 
239.233 
238.7 66 
287.557 
286.361 
285.226 
234.161 


299.638 

293.150 

291.277 

290.660 

290.205 

289.500 

289.260 

283.614 

287.606 

286.161 

285.202 

284.112 


*• Calculated temperature 
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(2) Conduction-Convection Model 

The errors inherent in the finite-difference approx- 
imations of .the jconducticn model are also present in the 
combined conduction-convection model. But the most- import- 
ant source of error is due- to hidden numerical dispersion 
terms,. These are terms that are introduced inadvertently 
into the original temperature equations (Equation 5) , when 
f inita=difference approximations that ignore second-order 
derivatives are used. In the backward- and the forward- 
difference approximations both the second-order space 
derivatives and the second-order time derivatives contri- 
bute to these dispersion error terms. In the central -diff- 
erence method, only the second-order time derivative 
introduces the numerical dispersion error terms. This kind 

of error cannot be avoided by using time-implicit finite* 

difference techniques. The numerical dispersion, errors 
are revealed in the central -difference method as follows:. 

If Equation (34i) is rewritten in terms of derivatives 
by using Equation (29) without ignoring the second-order 
time derivative, the resulting equation is: 

i^tT tt + DT/Dt * ^(T^ + T yy ) (42) 

By comparing Equation (42) and Equation (5), we find that 
an extra term totT^ has been introduced into the original 
Equation (5) by the finite-difference approximation that 
ignored second-order time derivative T^^. By differentiat- 
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i 

! . j 

| in G Equation~(5) with respect to t, ignoring third- and 

| ; higher order derivatives, and introducing the result into 

Equation (42) we get, on further rearrangement: 

DT/Dt = (^ L “^t)T n + U L -£v 2 At)T yy - uvAtT (43) 

Thus, the finite-difference (central-difference) approx- 
i imation has introduced the extra terms tT^,,., -;i-v 2 T^ t 

nnd “U^T X2 At. into the original Equation (5). These are the 
numerical -dispersion terms* The numerical-dispersion co- 
i efficients — £U At, -;?v At and -uvAt must be comnared against 

4 l to see how much error is introduced. They could introduce 
devastatingly significant errors. Errors due to dispersion 
' terms do not arise in the conduction model but in models 

involving the boundary-layer equations of energy and motion, 

* These errors can be reduced but not eliminated by including 

* i 

higher- order derivatives of time in the difference approx- 
imations or by reducing the magnitudes of the velocity 
components and/or time step. Table 3 gives the magnitude 
of the dispersion coefficients evaluated for our problem in 
which \ L = 9*306 x 10-* cm 2 /sec , u max « 13.92 x 10** 3 cm/sec, 
snd v 1Qax = 5.2 x 10"^ cm/sec. 

Note that, although it is desirable to decrease Ay 
and Ax and thereby increase maximum allowable velocities 
due to stability restrictions, the dispersion error terms 


4 i 


\ 



¥ 


*'■ *• 
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Table 3: 


kaomitijde, op Djrsppp^iorr^cp' 
to t;;v: t;: b : B al dTx¥j5iviT.: 


X^ = 9. 30 6 x.10”^ cm^/sec 


PFICIRMTS COMP ARED 
‘(x p of the LIQUID . 


u = 13.92 cm/sec; v = 5.22 cm/sec; . x = y = 0.127cm 


At 

&U 2 At 

§V 2 At 

(sec) 

(cm 2 /sec) 

(cm 2 / sec) 

6.0 

0.625X t 

0.038X t 

4..0 

0.^16X“ 

0.059Xr i 

2.0 

0.208X" 

0.029Xt 

1.0 

0.10^ 

0.015X7 1 

0.5 

0.052X" 

o.oo7x£ 


uvAt 

(cm 2 /sec) 

+0 • 468X T 
±0,312X7' 
±0.156Xt 
+0.078X" 
±0.039X£ 


increase with increasing velocity or decreasing Ax and Ay. 
We ran digital -computer programs for the conduction-convec* 
tion. model at different values of At for the. same Ax_=£y 
=0^127 cm that was used for the conduction model. We found 
that the dispersion terms did not. make any significant 
difference at At = 1.0 sec. At At _= 6.0 sec, the numerical 
dispersion terms made a substantial difference in the 


accuracy of the theoretical results. The time required to 
complete a computer program to acquire the equivalent of 
two hours of experimental time was prohibitive_for At ^ 1*0 
sec. Perhaps, a .faster computer would be very helpful. We 
conjecture that the numerical -dispersion coefficient -uvAt 
counteract* some of the effects of the coefficients -§u 2 ,At 
and -J-v At, so that, when the last two-mentioned are of the 
order of 10 percent of X L , the net error effects are much 
smaller than that. 
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It happened also, that in our system, the maximum 
temperature gradients in the liquid phase occur along the 
y-coordinate direction* that is, in the direction such 
that |vTyJ» luT^J. Thus, v had a lot more influence on the 
magnitude of the temperature profiles than u did. But a 
glance at Tatle 3 reveals that the error terms due to v 
alone in the dispersion terms are much smaller than those 
introduced by u. Thus, the numericaL-dispersion terms did _ 

not- have as much influence on the theoretical results as _ 

as they would if the maximum temperature gradients had been 
in a direction to cause JuT x |^> |vT | , 


expehikkiitai equipi;sxt a:jd pr ocedure 


Equipment 

The principal elements of the equipment were a tesJb 
cell, a 24-channel multipoint temperature recorder, and a 
refrigerator. The auxiliary elements were copper-constantan 
thermocouples, a power-driven liquid pump, methanol-, some 
pipes and tubings, and a test-cell stand that could be 
swivelled to various inclined positions from the horizontal 
plane. The test material, n-hexadecane , was liquid at the 
ambient temperatures for the experiments. 

Test C ell: The test cell (Pig. 1 and 5) had a con- 

stant cross-section of external dimensions 12.70 cm, and 
an overall height of 9.20 cm. It was composed of a copper 
cooling chamber soldered to one face of a 0.32-cm- thick 
copper plate (henceforth referred to as the cold .plate or 
the bottom plate) which was in turn bolted and glued to 
one end of a plexiglas frame. The frame was glued and 
bolted to a 1 .27-cm- thick plexiglas plate on the other end. 
Thus, the plexiglas frame, together with the plexiglas plate 
on one end and the copper plate on the other end formed a 
10.1 6 -cm-squar e , 3.81-cm-high cavity with 1 .27-cm-thick 
walls for containing the test material, normal hexadecane. 
The cooling chamber was constructed from 0.64-cm copper 
plates soldered together. 
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150 % of full size 
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The test coll carried twenty-four copper-constantan 


thermocouples located at certain measured locations. Table 
(4) details the thermocouple locations. The coordinates 
given in the table were located and fixed as shorn in Figure 
(1) and Figure (5). The positions could be off by, ± 0.16 cm 
because of the limitations in the measuring accuracy* 

Table THERMOCOUPLE LOCAT I ONS IN TEST CELL 

x * horizontal distance (parallel to the cold bottom 
plate), starting from left-corner of test-cell 
cavity. 

y = vertical distance from cold bottom plate— 

Each distance could be off by ± 0.16 cm 



No. « 

the number ass 

igned 

to a thermocouple 


No. 

X 

y v 

No. 

X 

y 

No. 

X 

y 


(cm) 

(cm) 


(cm) 

(cm) 


(cm) 

(cm) 

1 

2.54. 

0.00 

13 

2.03 

1.52 

3 

0*00 

2.54 

4 

8.13 

0.00 

20 

4.06 

1.52 

6 

10.16 

2.54 

14 

2.03 

1.02 

23 

6.10 

1.52 

8 

2.03 

3.05 

21 

4.06 

1.02 

16 

8.13 

1.52 

15 

4.06 

3.05 

24 

6.10 

1.02 

9 

2.03 

2.29 

18 

6.10 

3.05 

17 

8.13 

1.02 

19 

4.06 

2.29 

11 

8.13 

3.05 

2 

0.00 

1.27 

22 

6*10 

2.29 

7 

2.54 

3.81 

5 

10.16 

1.27 

12 

8.13 

2.29 

10 

7.62 

3.81 


Temperature Recorder: 

The 

recorder 

was a 

24-channel , 


multipoint recorder sold by ACCO Bristol , model 66A- 
24 PGC 570, that operated on 120 volts of 60-cycle alter- 
nating current. It recorded temperatures in degrees Fahren- 


4-1 

heit in the range -lQO °F to +200 °P ( 199.83 °K - 36 6.49 °K) . 
It had a chart drive speed of 2. 54. cm per minute and print 
speed of 2 seconds per point. The chart could be read with 
an accuracy of ±0.2? °K (0.5 °F) . 

Pump: The pump used to circulate the coolants (methanol) 

from the refrigerator to the test cell was a Chemical Rubber- 
Company ' J "No-Seal" centrifugal pump, Model ABIP005N#. 

It operated on 115-valts, 60-cycle alternating current only. 

It could attain 50 revolutions per second and pump from 
441.6 cc per second at a head of 30.5 cm to 262.8 cc per 
second at a head of 274.3 cm under normal atmospheric 
conditions* 

Refrigerator: The refrigerator for the coolant was a 

Bar Ray of Brooklyn, New York, Model 557T refrigerator that 
operated on a 60-cycle, 115-volt alternating current. It 
had a regulator that could be used to adjust the steady-state 
temperature to which the refrigerant was cooled. 

A schematic diagram of the assembled equipment is shown 
in Figure (6). A two-way tap-control valve permitted the 
circulation of the coolant in an auxiliary circuit -until 
temperature equilibration was achieved in the coolant. 

Experimental Procedure 

The test material, n-hexadecane, (liquid at normal room 
temperature) was introduced into the test-cell cavity through 
a filler-port. An excess of the material was allowed to 



FIGURE 6. BLOCK DIAGRAM OF ASSEMBLY OF MAIN EXPERIMENTAL EQUIPMENT. 
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collect in an expansion chamber connected, to a reservoir 
ol the -tes'C-material. In this-way, air pockets in the test 
cell were kept to a minimum as tno solidification progressed , 
by the introduction of new test material to assume the vac- 
ated volume. 

The temperature recorder was turned on at the same_time 
as the pump was turned on to circulate cold methanol to cool 
the bottom plate of the test cell. By use of the auxiliary 
circuit, a step change in the temperature of the test cell 
could be achieved. The experiments were performed with the 
cell inclined at the following angles from the horizontal 
plane: 0 °, 15°, 30 °, 4.5° and bu°. Studious attempt was 

made to approximate the same starting conditions for all 
runs — the same cooling rate, the same low temperature for 
the coolant and the same ambient temperature. The ambient 
temperature was_approximately controlled by setting the room 
thermos oat at a constant level for about 24 hours before 
an experimental run. Before each run the thermocouples 
were calibrated on the recorder by testing against room temp- 
erature as recorded oy a mercury thermometer. V/hen .all 
24 thermocouples recorded the same temperature within 
-0.27 °Iv, the experiment was begun. Room temperature was 
watched continually during the course of an experiment via 
a mercury thermometer. The temperature at the start of an 
experiment, as recorded by the 24 thermocouples, usually 
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did no ■' late from that recorded by tho mercury thermo- 
meter by ; ore than 1,1 °K, Plow of coolant was kept approx- 
imately constant from run to run by .opening all the valves 
to their fullest extent. The pump had only one speed setting. 
The low temperature of the coolant was regulated. by the 
11 cut in" temperature of the- refrigerator which was set- constant 
at about 272 °K. Methanol was used- as the coolant because 
it stayed liquid at this low temperature whereas water be- 
came ice. With such a setting, the coolant could beJbrought 
to a low temperature of about 270 °K. 

The experiments were usually terminated after about 
2 hours when more than half of "the test material had solid- 
ified. The temperatures, read from the charts, were then 
plotted against time. The bottom-plate temperatures were 
fitted into a polynomial function of time, comprising a ramp 
decay followed by a constant temperature. For a given set. 
of starting conditions* the experiments were repeated to 
test for reproducibility. 

The test material was practical n-hexadecane (n-C^H^) 
distributed by the Eastman Kodak Company^ 6 ) f or chemical- 
purposes. It had small impurities that influenced its solid- 
ification temperature. The solidification range was given 
by the manufacturer as 291.0 °I< - 289.2 °IC, However, we 
found that our sample froze at 290. b °K (17,5 °C). 


ai:d results 


SSS fe-£og-jl ho .Rcprod uclixni ty of Experl yg iental na<-.a 


Tv;o experimental rune under similar conditions of 
angle oi inclination of test cell, cut-in temperature of 
re rigor a tor, flow rate of coolant., chart speed and ambient 
temperatures as close as possible to each other we re te sjbe d 
for the reproducibility of the data. I t^was. estimated that 
experimental errors due to thermocouples, chart reading, 
and ambient temperature would cause as much as ±0.83_°K 
error in the experimental data. Thus, the- mean of the differ 
ences of- the temperatures recorded by the same thermocouple 
in both experimental runs were estimated with a 99 $ level 
of confidence. 


First of all the temperatures were plotted on the 
same temperature-versus-time graph. Then, 17 points were 
selected on the time axis and the temperatures of the two 
runs were read off • The difference of the temperatures 
was obtained for each of the seventeen points. Thus, a 
sample mean of the differences and a sample standard devia- 
tion were calculated. A confidence interval was calculated 
at the 99# level of confidence. If the confidence interval 
fell within ± 0.83 °K, the data were accepted as reproducible. 
These calculations were performed for each of the 24 thermo- 
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couples for a given pair of experimental runs under test. . 
Pig. 7 gives sample graphs for testing reproducibility of 

data for experimental Runs 9 and 10_, 

Let e^be the difference between temperatures recorded 
by the same . thermocouple between a pair of experimental 
runs under identical conditions at a point i in time and 
let e be the mean of a sample of size n of such differences. 
Then the true mean m Q of the differences, at 99 % level of 
confidence, lies in the interval 

e - ts/Vft < m e < e + ts/yfi inclusive (44) 

where s is the standard deviation of the sample of size n, 
and t is calculated from tne Student t-distribution at the — 
upper tail. 

For our system, ^ = 2.921, and n = 17^ 

n I o? - ( X e,) 2 

s » M Isi . W) 

n(n - 1) 

Thus, the confidence interval is 

e - 2.921s/Vl7< m e < e + 2.921s/Vl?. (46) 

This interval must lie within the limits of experimental 
error before a pair of runs may be accepted as being repro- 
ducible. The test for reproducibility of data was performed 


£ 


TEIIPERATURE, °K 






«, 4 
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for all 8 pairs of runs and for each thermocouple. Table 5 
shows a sample test result for experimental Runs 9 and 10, 


Table 5: REPRODUCIBILITY TEST BBTvJEEK RUN 9 AND RUN 10 

The decision, "accept" means that results, are repro- 
ducible, The confidence interval must lie within 
-0.83 °K for the results to be accepted as reproducible. 
Level of conf idence » 99$, 


Thermocouple 

Number 

1 

2 

l 

1 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 
17 

18 

19 

20 
21 
22 

2 
2 


e °K 

0.033 

0,065 

0,278 

0.033 

0.040 

0.251 

O.O69 

0.053 

0.125 

0.152 

0.180 

0.243 


s °K 

0.092 

0.251 

0.278 

0.092 

0.407 

0.334 

0.209 

0.458 

0.131 

0.316 

0.239 


Conf idence 0 
Interval K _ 


-0.033, 
- 0 . 112 , 
0.081, 
-0.033, 
-0.247, 
0.014, 
-0.079, 
-0.271 ,_ 
0.033, 
- 0 . 072 , 
0 . 001 , 


0.098 

0.243 

0.474 

0.098 

0.328 

0.488 

0.217 

0.377 

0.217 

0.376 

0.349 

0.420 


0.1.76 

0.154 

0.067, 

0.285 

0.229 

0.202 

0.086, 

0.372 

0.114 

0.167 

-0.004, 

0.232 . 

0_051 

0.117 

-0.032, 

0.1.34 

0.090 

0.343 

-0.153, 

0.333 

0.233 

0.244 

0.065, 

0.411 

0.066 

0.288 

-0.139, 

0.270 

0.040 

0.173 

-0.083, 

0.163 

0.033 

0.280 

0.082, 

0.231 

0.131 

0.210 

-0.018, 

0.279 

0.049 

0.293 

0.158, 

0.256 

0.107 

0.185 

0.024, 

0.238 

all 16 experimental 

runs in 

pairs a 


to reproducible pairs. 


Decision 

accept 

accept 

accept 

accept 

accept 

accept 

accept 

accept 

accept 

accept 

accept 

accept 

accept 

accept. 

accept 

accept 

accept 

accept. 

accept 

accept 

accept 

accept 

accept 

accept 


JfesBBa 
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Table o; 


EXPERIMENTAL PUNS LISTED 


IN REPRODUCIBLE PAIRS 


First values of temperature correspond to experimen- 
tal runs shovrn first in the first column* 


Runs 

<x° 

1; 

2 

0O 

3; 

4 

00 - 

5; 

6 

15° 

7; 

8 

600 

9; 

10 

300 

11; 

12 

45° 

13; 

14 

45° 

15; 

16 

60° 


T a , °K- Cold plate steady 

temperature s . °K ' 


300.7 

300 . 7 

270.4 

2.7u.^ 

300.1 

300.4 

269.8 

269.8 

298.7 

298.7 

270.4 

270.4 

298.7 

298.7 

270.4, 

270.4 

301 . 2 : 

1 301.5 

269.8: 

269.8 

300 . 1 ; 

; 300.1 

269.8; 

269.8 

301.5; 

; 301.5 

269.8; 

269.8 

301.5; 

; 301.5 

270.4; 

270,4 


Presentation of Experimental and Theoretical Data 

The experimental data_for all the 16 experimental 
runs are presented in Appendix A. Figures (8) to (18) show 
comparisons of experimental and theoretical data. The 
experimental runs were performed under conditions listed 
in_ Table 6. Figures (8) and (9) were obtained for experi- 
mental runs with the cell sitting on a horizontal plane. 
Remarkably good agreement was obtained—between the experi- 
mental data and theoretical conduction data, thus indicating 
that when the cell was sitting on a horizontal plane, heat 
transfer was by conduction with negligible convection. 
Obviously, a convection model could not predict the results 
of Figures (8) and (9)» The data have been presented accord- 
ing to distance y from the cooled bottom plate. The data 
from the experimental runs at angle 0° (horizontal) indi- 
cated that temperatures were functions of time and distance, 
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y*_££pni the cooled plate, out not of horizontal distance x, 
as 'would be expected from a conduction model with our given 
boundary conditions. 

Figures (10) to (18) present results obtained with 
the test cell inclined at different angles from the horizontal, 
starting with a 15° angle. A study of the experimental points 
on these graphs reveals the presence and importance- of con- 
vective heat transfer. Now, the temperature profiles .become 
functions of horizontal distance x. Here, again, as in 
Figures (8) and (9 )» the temperature graphs have been arranged 
according to vertical distance y from the cold plate. At_a 
given vertical distance y, the maximum and minimum tempera- 
tures along a horizontal direction, x, have been plotted on 
the same graph. Other thermocouples with in-between temper- 
atures have been plotted singly. Thus, the separation- of 
temperatures caused by convection along a horizontal direction 
was revealed in each case. 

The general effect of convective heat transfer would 
appear to be increased heat-transfer rate, resulting in lower 
temperatures throughout the cell and faster freezing rate 
than would be obtained with conduction heat transfer only* 

Also, convection changed the shape of the interface from 
being flat to a ourved shape (Fig. 13). More freezing occurred 
near the pivot x = 0 of angle of inclination than near x * W. 

The effect of increasing angle was to increase the 
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cooling-down of that part of the cell to the right of the 
plane x = W/2. The greater the angle of—inclination, the 
more the cell_is cooled down below pure conduction tempera- 
ture- profile. It would seem that after certain temperature 
gradients were achieved, the cooling rate was greatly re- 
duced. The magnitude of maximum deviation of temperature 
along a horizontal direction at a given height, y, from the _ 
cold plate was observed experimentally to depend more-on the 
temperature difference between initial temperature T and 
solidification temperature, T f . Rough estimates of this 
difference were found to be approximately £(T - Tu>) Irre- 

spective of the angle of inclination .—However, the rate of 
attainment of such maximum differences differed with angles 
and could not be estimated. 

Figures (10) ...to (18) show trends in temperature pro- 

files as solidification with convection progressed. The 
theoretical curves could not be brought to match closer 
with the experimental curves, because stability limitations 
from the finite difference approximations of the tempera ture 
equations would not allow us to use the maximum possible 
velocities from the velocity profiles. 

If velocities could be increased without violating 
the stability criteria of the finite difference solutions, 
then it is felt that it would be possible to reach the proper 
maximum velocity when the theoretical and experimental curves 
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would be ~.ry close to each other. At such a time one could 
calculate empirically the maximum velocities for the system 
as functions of angle a. VJith the true maximum velocities 
known, one could then find the proper value_of (T - T) 
involved in the gravity approximations in the velocity 
profiles. 

Rough estimates of fractions of maximum velocities 
that we actually used were obtained as follows. At the 
start of convection, a critical Rayleigh number was assumed 
based on values available in the literatures Prom the 
critical Rayleigh number, AT at the start of convection vras 
computed and used in Equation (28) to calculate the maximum 
velocity component u# max ( > v rnax ) . Then, ratios of the 
actual u^. used and this computed u** max were .obtained for 
different angles. .The ratios were then the fractions of 
initial maximum velocity (or in other words, gravity levels) 
that we could attain at the very start of convection. Equa- 
tion (4?) gives the relationship for the velocity ratios or. 
gravity level Q: 

Q * 0.95 (2X L /Ax)/B or * 45.6h ( /R ax sin a — _(4?) 

where R acr ® critical Rayleigh number 

* - SGh 0 3 JT 0or / ( ) 

and h Q « height of liquid in test cell at the start of 
convection . 
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A half of the height of the liquid phase at the start of 
convection was used because the temperature gradient was 
still confined to y< h/2. The factor O.95 appeared in 
Equation (47) because we used 95/' of the maximum velocity 

allowed by stability requirements of the finite difference 
solution. 

Table 7 shows gravity levels (fractions of maximum 
velocity obtained) against magnitudes of angle of inclina- 
tion oc . It shows that, as the angle of inclination in- 
creased, our stability- allowable maximum velocity became 
a smaller fraction of the maximum velocity predicted by 
Equation (28) at the start of_convecti;,n. As solidification 
progressed, the velocity in the system could approach the 
actual-maximum velocity that we used and could also- become 

smaller than that. It must be emphasized that the gravity 

levels in Table 7 correspond to the start of convection only. 


7 • GRAVITY LEVELS AT START OP CONVECTION, (q ) 

Jx = 0.12? cm, R = 3,500 

actuall Y used — l*9\j^/ ax — 13*92 X 10**^ cm/seo 
X L * 9.306 X 10* 4 cmVsec 


go 

15 ® 

30° 

45 ° 

60 ° 


Gravity Level. Q. 

0.68 

0.39 

0.28 

0.23 


PERATURE, Ok TEMPERATURE 


^UURE 8. EXPERIMENTAL AMD THEORETICAL DATA FOR RUN_2: 
( * ~ 0° ; T a « 300.66°K 


(a) Temperature of cold plate (y=0), 

x Thermoe'ple Expt. Conduction 
(cm) No. (theory) 
2.54 


8.13 


TIME (100 SEC) 


(b) Temperature at yal.02 cm . 

x(cra) Thermoc. Sxpt* Co^gJ^on 




if i 


4 3 

E (100 SEC) 






TEMPERATURE, °K TEMPERATURE 


run 2 : a » o° 


300.66 °k 


(c) Tempera ture ab y •* 1,2? cm 


x Thermoc. 2::pb. Conduction 

0*00 2 O (the “’ y) 

10.16 5 o 


TIME (100 SEC) 


(d) Temperature at y = 1.52 cm 

x- Thermoc. Expt. Conduction” 
(cm) (theory) 

2.03 13 £> 

4.06 20 0 

d 6.10 23 A - 

\n 8 ‘ 13 16 O ZI 


TIME (100 SEC) 

— 1 1 j i 

(e) Temperature at y = 2.29 era 

x^(cm) Thermoc. Expt. Conduct (theo' 

19 8 

22 ^ 

12 © 


irAOKi] 


TIME (100 SEC) 












IMPERATIVE 


:oup.3 8 (contd). RUM 2 


300.66 °K 


f) Temperature at y 
x ( c i.i ) Thermoc. 
0.00 3 

10.16 6 


= 2.54 cm 

Expt. Conduct (theory) 


TIKE -(100 SEC) 


(g) Temperature at y = 3.05 cm 

x(cm) Thermoc. Expt. Conduct (theory) 
2.03 8 © 

4.06 15 Q — 

6.10 18 A 

8.13 11 0 ~ 



TIKE (100 SEC) 


(h) Temperature at y = 3. 81 cm 

x (cm) Thermoc. Expt. Conduct (theory! 
2.54 7 © - 
7.62 10 Q ~ 



TIME (100 SEC) 









HATURE. °K TEMPERATURE 








I(. JUS 9 ( contd) . JIU.T 3: a 


o 


300.11 °K 


(c) Temperature at y » 1.2? cm 

"(cm) Thermoc. Expt. Conduct, (theory), 

0.00 20 

10.16 5 0 


(dl_Temperature at y = 1.52 cm 

x(cm) Thermoc. Expt. Conduct^ theory 
2.03 13 O 

4.06 20 □ 

6.10 23 A 

1 8.13 16 q 


(e) Temperature at y = 2.29 cm ' " 

x ( cm ) Thermoc . E^pt . Conduc t ( theory 


5*0 


2 


3 






9 (oontd) . Eli: I 3: n 


'•■onporatm’c at y - 2 , 54 cm 

x(cm) Thermoc . Expt . Conduct (theory) 


TIME (100 SEC 


(g) Temperature at y ■ 3.05 cm 

x(cm) Thermoc. Expt, Conduct (theory) 
2.03 80 

4.06 15 0 

6.10 13 & 

8.13 11 o 


TIME (100 SEC 


(h) Temperature at y = 3,81 cm 

x(cm) Thermoc. Expt. Conduct (theory) 
2.54 70 

7 Jl 2 . 10 0 muzz 




& © 


TIKE (100 SEC) 












TEMPERATURE 


FIGURE 10 . EXPERIMENTAL AMD THEORETICAL DATA FOR RUN 6 

<* “ 15 ° ; T a = 298.72 _ 2 K 


(a) Temperature of cold plate (y = 0 ) 

x(crn) Thermoc. Sxpt. Conduct (theory) 


(b) Temperature at y =1.02 cm for x=2.03cm 
and x= 8.13 cm. 

x(cm) Thermoc. Expt. Convection (theo) 
2.03 lty O — » •* — • — • 


8.13 


<£v.o 


17 ._0 

— Conduct(theo) for all pts. 




©o 




0 ° OS' 


IHE (100 SEC) 









PERATURE, °K TEMPERATURE, °K TEMPERATURE 


FIGURE 10 (coafcd) 


298.72 °K 


(oJ_JCfi.mporature at y = 1.02 cm for x=4.06c 

Thornoc. Expt. Convection( theo. ) 

21 0 — — — 

n Conduction (theo) 


TIME (100 SEC) 


(d) Temperature at y =.1.02 cm for 
x = 6.10 cm 

Thermoc. Expt .Convection( theo- ) 

) 24 0 — — — 

v 0 

— Conduct ion (theo.) 




TIME (100 SEC) 


(e) Temperature at y = 1.27 cm 

x(.cm) Thermoc. Expt. Convection (theo) 

0.00 2 O — • — . — 

10.16 5 0 


Conduction (theo.) 
for all pts. 








TIME (100 SEC) 













PERATURE, °K TEMPERATURE, °K TEMPERATURE 


FIGURE 10 (co'Atd) . RUN 6: a = 15° ; T a = 298,72 °K 

305 p — — 1 1 

(f) Temperature at y = 1.52 cm for 
x = 2,03 cin and x = 8.1 3 cm. 

oao x(cm) Thermoc. Expt. Convect. (theo. ) 

JUU - 2.03 13 © — . — 

8 « 1 3 16 O — 


j&N?© c< 

®3k§La o : 


Conduction ( theo. ) 
for all pts. 


r<m 

ss 


12 24 36 So 

TIME (100 SEC) 

' ' .. 

(g) Temperature at y = I.52 cm for 
x = 4.06 cm 

^ Thermoc. Expt. Convection (theo.) 


Conduction ( theo . ) " 

0 G / 


TIME (100 SEC) 


(h) Temperature at y = 1.52 cm for 

x = 6.10 cm. 

Thermoc Exjj>t. Convection (theo.) " 
Conduction (theo.) 

— ii! A ^ 



TIME (100 SEC) 







TEMPERATURE , °K TEMPERATURE 


IGURE 10 (conUl). RUN 6: a = 15° ; T a = 298.72 °K 

| - -* — 1 1 L 1 ■ 

M 305 - (i) Temperature at y = 2.29 cm for 
6 2: = 2.03 cm and x = 8.13 cm. 


_ G-v* 


x(cm) Thermoc. Expt. Convect (theo.) 

2.03 9 O — . 

8.13 12 0 — — _ 


Conduct (theo) 

(SUv /-x for a11 pvs. - 


m 


TIME (100 SEC) 

1 j ■ - 

(j) Temperature at y = 2.29 cm for x=4.06 cm 

Thermoc. Egpt. Convection (theory) 
Conduction (theory) 


. 00 © 


0© 


24 36 48 

TIME (100 SEC) 



TIME (100 SEC) 





IGtTRK 10 (contd) . RUN 6: a = 15° ; T_= 298.72 °K 


305 r (1) Temperature at y~2»5^ cm. 


x(crn) Thermoc. Expt. Convect (theo.) 

0.00 3 O — • — • — 

10.16 6 0 — - — — 


r °o- 

I 295 - °\ 


-Conduction (theo.) 
for all pt§ 


GOO O' 


2 h 36 W 

TII1E (100 SEC) 


(m) Temperature at y = 3.05 cm for 
x. = 2.03 cm and x = 8.I.3 cm. 


x(cm) Thermoc. Exot. Convect (theo.) 
2.03 3 ©‘ — . — . — 


8.13 11 0 — 


Conduction (theory) 
for. all pts. 




TIME (100 SEC) 


(n) Temperature at y=3.05 cm for x=4.06 cm. 


Thermoc. Expt. Convection (theory) 

15 Q — 


Conduction (theory) 




''©~00 . 


1 Q- -O 0 ». Q_ 0)_ Q 


(100 SEC) 




■rm I 








TEMPERATURE , °K TEMPERATURE 


FIGURE 10 (contd) . RUN 6: a 


293.72 °K 


_ Temperature at y=3.05cii for x=6.10 cm 

Thermoc. Eg)t. Convection (theory) 

Vj. Conduction (theory) 

- ® 

~ 1 p- _ 

12 24 3 6 48 60 \ 

TIME (100 SEC) 

« . 

(p) Temperature at y=3.81 cm 

x(cm) Thermoc. Expt. Convect (theo.) 

2.54 7 O . — 

7.62 10 0 — — — 

Conduction (theo.)_ for 


24 36 48 

TIME (100 SEC) 


FIGURE 11. IIIITIAL STREAMLINES IN TEST CELL (THEORETICAL) 

3. Sir -t- 


o.oerT 


0=0.0044 cra2/sec 


2.54 5.08 7.62 10.16 

DISTANCE PARALLEL TO COLD PLATE (x cm) 






SOLID HSIC-HT (Y cm) ^ SOLID HEIGHT 


FI GLEE 12. HEIGHT OF SOLID PI: AS 


x(cn) Expl*,. Convection (theo.) 

2.03 0 — — — 

4.06 ® — * — • — — 

0.13 A — + — t — 


Vik 1 


C on&ucfc ion ( theo. ) 
for .all x 1 s. 


TIME (100 SEC) 


FIGURE 1*3. INTERFACE SHAPE AT t = 7200 SEC (THEORETICAL) 










TEMPERATURE 


FIGURE 14 


EMPERILEHTAL AMD THEORETICAL DATA FOR RU1J 7; 
a = 60° ; T a = 293.72 °K 


— - ■ — I 

^ (a) Temperature at 

- 1 ‘ 

y = 0 (cold plate) 

_ x(cn) Thermoc. 

Conduction _( theo ) Axpt • 

2.54 1 

0 " 

8.13 4 . 

0 

1 ) 

tm 

-%-0 Q 0 Q 
^ ■ , — 

6 0 0 0 8 0 0 


t- 3 6 48 

TIME (100 SEC) 


(b) Temperature at y = l o 02 cm for 
x = 2.03 cm and x = 8*13 cm. 

x(cn) Tliermoc. Expt. Convect (theo.-) 

3 °°r 2.03 14 0 — * 

0 8.13 17 0 

Conduction (theo.) 

295 for all pts. 


TIME (100 SEC) 




•IPERATURE, K TEMPER ATUBE , K TEMPERATURE 


[»T»] 


FIGURE- 14 (contd) . RUM 7 


293.72 °K 


(c) Temperature at y=1.02 era for x=4.06 cc 

Thcrraoc. Expt. Convection (theoJ 

21 O — — — 

— Conduction (thco.) 


O <5 


tii:e (100 sec) 
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FIGURE 15 (contd). RUM 9; a = 30° • T Q » 301.22 °K 
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GURE 15 (contd). RUN 9 : a = 30 ; T a = 301.22 °K 
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FIGURE 16. EXPERIMENTAL AITD THEORETICAL DATA FOR RUN 11 
a = 45° ; Ta = 300.11 °K 
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FIGURE 16 (contd) 
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FIGURE 17. EXPERIMENTAL AND THEORETICAL DATA FOR RUN 14* 
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FIGURE I? (contd)* RUN 14: a » 45° ; T a « 301.4.9 °K 
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MGURE I? (contd). RUi! 14 : a = 45 ° ; T a «* 301.49 °K 
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FIGURE 17 (contd). RUN 14: a 
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FIGURE 18. EXPERIMENTAL AMD THEORETICAL DATA FOR RUN 15: 
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FIGURE 18 (contd). RUN 15: a=60° ; T a = 301.49 °K 
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FIGURE 18 (contd) . RUN 15; a 
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FIGURE 18 (contd). RUN 15: a = 60° ; T a = 301.49 °K 
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FIGURE 10 (contd) 
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FIGURE 13 (contd) * RUN 15: a ■ 60° ; = 301.49 °K 
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CONCLUSIONS AND REGOliKENDATIOUS 


The following conclusions were drawn from this study: 

1. When the test cell was sitting on the horizontal 
plane with angle a = -Q_°, the heat transfer mode was by 
conduction, not convection. 

2. Good agreement was obtained between the pure- 
conduction model and the experimental rims made at a = 0°. 

3. When conduction was the_heat transfer mode, the 
shape of the solidification interface was planar., not _curved. 
For all horizontal x- positions the solid-phase height was 
the same . 

4. The temperature during conduction solidification 
was strongly a function of time and distance y from the 
cooled bottom plate. It had no strong- dependence on the 
horizontal x-direction when the containing side walls acted 
as insulators. It was found that one-dimensional transient 
conduction model in the y-direction predicted the phase 
change as well as did transient two-dimensional model that, 
included the x-direction. This finding supports the find- 
ings in a previous transient one-dimensional conduction 
phase change study made by these investigators (Reference 13). 

5* In solidification by conduction, in which the side 
walls approximate insulators, the extra minor accuracy ob- 
tained by using a transient two-dimensional model was not 


97 


93 


worbh the large differential in cost in computer time from 
a one-dimensional transient solution. A uni-dimensional 
transient solution would do as well and still save computer 
time. 

6. The phase-change_process was the controlling 

process. Whether the heat transfer mode was by conduction 
or convection, the actual values of phase-change enthalpy 
change, ..solidification temperature, density, specific heat 
and thermal conductivity at the interface seemed to control 
the accuracy of prediction of experimental data by theoretical 
models. More accurate values must be found for these proper- 
ties at the interface where some .Of them experience disconti- 
nuities and sudden jumps. In the theoretical calculations, 
the interface was found to act as a pseudo-insulator between 
the two phases* present. 

7. With the container of the test material inclined 
at an angle, we found that convective effects became impor- 
tant. 

8. Convection increased the heat transfer rate and 
caused a faster overall cooling of the test cell and its 
contents. Two indications of the magnitude of convective 
effects in the cell were how much below the conduction temp- 
erature the cell had been cooled at any given time, and the 
shape of the solidification interface* With convection 
present, and cooling occurring from below, the interface 
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became shaped slightly like an elongated S-curve with the 
up branch near the vertical wall passing through the origin 
of the angle of pivot, the flat (or inflexion) part near the 
half-way mark between the two side walls. The shortest 
height of solid occurred near the far side wall . 

9. We could not get a complete solution of ‘the full 
velocity equation by finite difference methods because of 
stability restrictions imposed by the solution procedure. 
Approximate velocity profiles were substituted into the 
temperature equation. The profiles involved a single cir- 
culation flow symmetric in the vertical y-direction but 
off-centered in the horizontal x-direction. 

10. Convection caused the temperature profile to 
depend strongly on both horizontal and vertical directions 
(unlike .conduction) even when the side walls were insulated. 

The maximum deviation in temperature between two points in 
the same horizontal plane at a given vertical distance fr o m 
the cold plate seemed to depend more strongly on the temper- 
ature difference between the ambient and the freezing interface 
than on the angle of inclination* The angle of inclination, 
however, had significant effect on the overall heat transfer 
process. 

11. Solid-solid phase transition may take place during 
solidification and the physical properties of the system 
should be modified to account for this. 
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12 » The theoretical calculations used in this study 
to investigate natural convection are first approximations 
and clearly delineate the trend and significance of convec- 
tion, With the availability of a much faster computer with 
a much larger memory bank than the computer, available to us* 
one may finally use small enough time steps and space incre- 
ments that would yield velocities higher or equal to those 
predicted analytically and still stay in the stable regime 
of finite difference solutions. In such a case, one_could 
estimate more accurately the shape and magnitude, of the 

velocities and come up with better gravity level approxi- 

roations in the velocity equations, 

13. Many theoretical methods used in the literature 
could not stand up to the test in predicting actual experi- 
mental data. These theoretical soluticns_were obtained 
under certain conditions and with physical properties that 
could not.be duplicated easily in the laboratory, 

14. Perhaps, a better approach to phase-change 
thermal control would be to use polymers that undergo solid- 
solid transition with high enthalpy change. Thus, the phases 
change material could stay solid or amorphous within the 
range of its temperature of operation. Thus, better packag- 
ing would result. Materials that melt or solidify are more 
difficult to package in the liquid phase and more difficult 
to predict as to performance because of convective and other 
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effects. 

15. Another approach to the thermal control problem 
with materials that so from liquids to solids- and vice versa 
would be to come up with empirical effective thermal proper- 
ties that incorporate convective and conduction heat transfer 
with phase change. For instance, such a procedure could 
take into consideration entrapped air pockets and voids in 
the solid phase, and heat transfer coefficients due to con- 
vection in the liquid phase. 

We conclude that, with cooling occurring from below, 
gravity-induced convection is important in the solidification 
of n-hexadecane when the cooling is not oriented perpendicular 
to the direction of gravity. Heat transfer rate.is, in 
general, increased. 


NOMENCLATURE 


Text Computer 
B A1 

c 

C p CPL , CPP , CPS 

°pS cp • 

d,d ROL,ROP-,ROS 

^ D(I) 

e i 


f(£) T(I,1) 

S 

G 1" G 11 G1 “ G11 
h AHO, AHI 

h p AHP 
H f HP 


Constant (Eq.24c) — cm/sec 

Constant (Eq.15) — sec-°K/cm 

Specific heat — cal/(gm-°K) 

Effective specific hest of- solid- 
phase (Eq.38) — cal/(gm-°K) 

Density — gm/cm ^ 

Coefficients in polynomial -fit, 
f(t) = D^Dgt+D^t^. ,+D i t i “ 1 of cold 
plate temperature — °K/(sec*“‘ ! ') 

Difference in temperatures recorded 
by the same thermocouple at the same 
time in a pair of experimental runs 
made under- identical conditions - °K 

Sample mean of e^ in a sample of 

size n — °K 

Temperature of cold plate — °K 

Acceleration due to normal gravity 
cm/sec^ 

Coefficients defined in Eq.30 and . 
Eq.31 

Height of the liquid phase of the 
test material (y direction) — cm 

Thickness of walls of test cell — cm 

Latent heat of solidification — cal/gm 
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Text Computer 
k AK,PK,SK 

K 

m 

e 

M M 

N N 

NS(i) NS (I) 


PtP 

q 

Q 

R a 

s 

t TIME,TSE 

At DELT 

T T, TO 

T a TA 

T f TP 

u UX 

V VY 


Thermal conductivity — watt/(cm-°K) , 

cal/(cm-sec-°K) 

Degree of polynomial fi±_f (t) of cold- 
plate temperature 

True .mean of e^ — °K 

Number of finite-difference, space nodes 
in the coordinate x direction 

Number of flutter-difference space nodes 
in the _ coordinate y direction 

Number of finite-difference space nodes 
in the solid phase in the y direction 

Pressure 

Characteristic number (Eq.l^) — cm“^ 

Fraction of normal gravity (Eq.47) 

Rayleigh number (Eq,^7) 

Sample standard deviation of e^ in 
sample of size n (Eq«^5) — °K 

Time — sec 

Finite-difference time increment — sec 

Temperature — - °K 

Ambient (room and initial) tempera- 
ture . — °K 

Equilibrium solidification tempera- 
ture of test material — °K 

Coordinate-x component of velocity - 
cm/sec 

Coordinate-y component of velocity - 
cm/cec 


Text Computer 


W AL 

x XA,XP 

Ax DELX 

X(x) 

y X 

Ay DELY 

a 

3 

Y BB 

X ALA, ASA 

v 

CO 

9 PHI 

it PI 


Width of cell cavity (x direction). - 
cm 

Abscissae in a Cartesian coordinate 
system — cm 

Finite-difference space increment in 
the coordinate-x. direction — cm 

Function defined by Eq.18 to Eq.22 - 
cm/sec_ 

Ordinate in a Cartesian coordinate 
system — cm 

Finite-difference space increment in 
the coordinate-y direction — - cm 

Angle of inclination of the test cell 
measured anticlockwise from the_hor- 
isontal plane — degrees 

Coefficient, in equation of state for 
density: d L _= d L(J + 0T — gm/(cm^-°K) 

Constant (Eq*23) ~ cm*’'*' 

Thermal diffusivity — cm 2 /sec 

Kinematic viscosity — cm 2 /seo 

Vorticity (Eq.10) — sec" 1 

Stream function (Eq.9) — cm 2 /sec 

Constant_= 3*14159266 


Subscripts Subindices 


1 


I Identifying number for finite-diff- 

erence node in the x direction 

Identifying number for finite- 
difference node in the y direction 


J 


S ubscript s 


L 

P 

S 


Liquid phase of test material 

Plexiglas, material of walls of test 
cell 

Solid phase of test material 


Superscript 

• Implies "evaluated at a new time " 

(Eq.32 to Eq.3^) 
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APPENDIX A 


The experimental data for the report is given in 
the following reference: 

Ukanwa, A. 0. , "Thermal Modeling of Phase Change Solid- 
ification in Thermal Control Devices Including Natural 
Convection Effects," Thesis No. T1422, Colorado School 
of Mines, Golden, Colorado, 1971 


APPENDIX B 


PORTRAIT IV Computer- Program ll C0MDET.F4 11 - for solving 
either the conduction model or the combined conduction- 
convection model of solidification heat transfer. 

This.. program was actually run on a P.D.P-10 digital 
computer via a time-sharing teletype. 

Instructions For Th e Use of "C0NDET.F4 11 • 

A. For program execution: 

6= INPUT DEVICE 7=0UTPUT DEVICE 

B. INPUT DEVICE (OR INPUT FILE) 

The cold-plate (bottom-plate) temperature is assumed 
to be a polynomial in time (t) for 0 < t <. TB, and constant 
thereafter at a temperature TC. Follow the following steps 
(in order) in reading in the input data: 

1. Consult the section on nomenclature in this text 
for terms used in the computer program. Then read in: 

2 * Line 1 ( or card 1 ) : Ambient temperature ( or initial 
uniform temperature TA) ; equilibrium temperature of solid- 
ification (TP); final constant cold-plate temperature (TC); 
finite-difference time increment (DELT); space increment in 
the x-direction (DELX); space increment in the y-direction 
(PELY); all in that order. Use free-floating point format of 
FORMAT 2 , 

3. Line 2 (or card 2): Time at which cold-plate temp- 
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erature becomes constant (TB) ; time at which program exe- 
cution should be terminated (TAU).Use free-floating point 
format (FORMAT 2)_.’ 

4. Line 3 (or card 3) : Degree of polynomial fit of 
cold-plate temperature (K) ; number of finite-difference 
nodes, in the x direction (H), in the y direction (N); an 
integer (NT) that is negative or zero for conduction-model 
solution but positive non-zero for the combined conduction- 
convection model solution. Use free-integer format (FORMAT 3) 

5. Line 4 (or card 4): Coefficients (D(I)) of poly- 
nomial fit to cold-plate temperature. (See the section on 
nomenclature for the definition of D(I) or-D^). Use free- 
floating point format (FORMAT 2) . 

6. Line 5 (or card 5) L_ Liquid-phase thermal conducti- 
vity_XAK)j solid-phase thermal conductivity (SK) ; liquid- 
phase density (ROL) ; solid-phase density (ROS) ; liquid-phase 
specific heat (CPL); solid-phase specific heat (CPS); latent 
heat of solidification (HP), Use free-floating point format 
(FORMAT 2). 

7» Line 6 (or card 6): Properties of material of 
the wall of the test cell: thermal conductivity (PK); density 
(ROP); specific heat (CPP). Use free-floating point format 
(FORMAT 2). 

8. Line 7 (or card 7) • Initial height of test-cell 
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cavity (AKO) ; width of test-cell cavity (AL) ; thickness 
of test-cell ^alls (AHP) . Use free-floating point format 
(FORMAT 2). 

A sample input file looks like this: 


301.49.290.6.269.83.1.0. 0.254.0.127 

362.0. 720e.0L 

1,81.31,3 , 

301. 49 ,- 0 . -0874 

0.00036,0.00037,0.755,0.833,0.506,0.501,56.67 

0.000496,1.155,0.35 

3.81,10.16,1.27 

Note that the listed input file is for the combined con- 
duction-convection model, since HT=3 is positive and non-zero. 
Also, according to this file, the cold-plate temperature is 
linear with time (K=l) until t=TB=362.0 when it stays constant 
at TO 269 . 83 . 

C. OUTPUT FILE 

This program _will print out 

1. The time at which a node in the y-direction solid- 
ifies, the node itself and its horizontal location. 

2. The temperature profile and the stream function 
profile of the entire test cell every 60 time steps. It starts 
with the cold plate (J=l) , prints out all horizontal (I) values 
and then goes to J=2, and so on. 

3. Adequate headings are printed to identify what is 
being printed out. 

The print-out may be modified from the general format 
given In this program to any desirable format by changing 
lines 64 to 78 of C0ND5T.F4 as needed. 
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APP?-?PIX_C_:_ OTHER THEORETICAL CQESTnERATTnM.q 


If one attempts to solve the vorticity equation 

(Eq.ll) or the velocity equations (Eq.8) by central-differ- 

ence methods,, one finds.. that the following stability re- 
quirements must, be satisfied; 

* 

1 - 2^t(l/(Ax) 2 + 1/ (Ay) 2 ) > 0 (C-l) 

l u (i» J)l < (C-2) 

I v(i, d)| < 2^y (c-3) 

When the vorticity or the .velocity equations are solved 
conjointly with the temperature equation (Eq.5), the con- 
ditions of Equation (C-l) to Equation (C-3) must be satis- 
fied conjointly with the conditions of Equation (39) and 

Equations (41). The most restrictive conditions are used 

to evaluate, maximum time increment (/it) and velocities allowable 
by stability requirements. 

For our test material, v >\ L ; hence, Equation .(C-l) 

rather than Equation (39) would determine At , but the 

rustic * 

allowable maximum velocities would be determined by Equation 
(41) rather than by Equation (C-2) or Equation (C-3). 

At the first time step, every term except the gravity 
term in Equation (11) or Equations (8J would be zero. 

Thus, from Equations (8), the velocities obtained at the 
first time step would be: 


u'(i,j) * - _ 3 . g6t ( T(i, j) -T ) sin a 


( C— 4) 
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v'(i,j) = - ( T(iJ)-T ) cos a 

\ 

The values given by Equations (C-4). and (a-5) must satisfy 
the condi tions_C.or stability given by. Equations (4-1) , (C-2) , 
and ( C— 3l_*. But, with experimentally-observed values, the 
temperature difference, T(i,j)-T, is such that the stability 
requirements of Eq.uations-441) , (C-2), and (C-3) are violated 
even for At_= 0.05 sec and Ax = O.O635 cm. It is very hard to 
complete stable computer solutions for values smaller than 
these. Yet, in order to get stable solutions of the velocity 
and the temperature equations at the very first time step, 
one needs a time-step size much smaller than 0.05 aec._Even 
with such small time steps, it would get more difficult to 
satisfy the stability requirements as the velocities would 
increase with subsequent time steps. This was the problem 
that forced us to use approximate velocity profiles for 
our study in stead of solving the complete velocity equa- 
tions by finite-difference methods* 


